Rydberg crystallization detection by statistical means 
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We investigate an ensemble of atoms which can be excited into a Rydberg state. Using a disor- 
dered quantum Ising model, we perform a numerical simulation of the experimental procedure and 
calculate the probability distribution function P{M) to create a certain number of Rydberg atoms 
M, as well as their pair correlation function. Using the latter, we identify the critical interaction 
strength above which the system undergoes a phase transition to a Rydberg crystal. We then show 
that this phase transition can be detected using P{M) alone. 

PACS numbers: 32.80.Rm, 75.10.Pq, 64.70.Tg, 34.20Cf 



I. INTRODUCTION 

Recent experimental progress in producing and con- 
trolling highly excited atomic and molecular aggregates 
have triggered a number of experimental and theoretical 
works on interacting Rydberg systems ii"— One of the 
most prominent phenomena observable in such systems is 
the dipole blockadejiiii^ which is a consequence of effec- 
tive interactions between atoms in Rydberg states with 
principal quantum numbers 71 ~ 30 — 80 J^ The physical 
reason for the blockade can be summarized as follows: 
the large dipole moment of a Rydberg atom induces size- 
able energy level shifts in the atoms in its vicinity. As a 
result, atoms within a certain blockade radius Rb cannot 
be excited, even though they arc subjected to the same 
electromagnetic field as the proper Rydberg atoms. 

Typical experiments are conducted on atom ensem- 
bles at ultralow temperatures and probe these systems on 
timescales during which almost no particle movement is 
possible. Therefore, the cloud of atoms can be considered 
as 'frozen' in a more or less disordered constellationiiSS 
After the required fine-tuned electromagnetic fields are 
switched on, a number of atoms undergo a transition 
into the highly excited Rydberg states. If Rb is larger 
than the average interatomic distance, only a fraction of 
the atoms can be excited, while the rest remains in the 
ground state due to the blockade effect. 

The spatial arrangement of the Rydberg atoms within 
the cloud has very interesting features. In Ref. [21[, the 
concept of a Rydberg crystal was put forward. The block- 
ade region formed around an excited atom can be mod- 
eled as an effective repulsive interaction between the Ry- 
dberg atoms. It might be responsible for an emergent 
long-range order of the Wigner crystal type^. However, 
its detection is extremely difficult. The primary method 
for detecting long-range order is spectroscopy, which is 
difficult to reliably realize in experimentsi^ Other ex- 
periments benefit from the low ionization energy of the 
Rydberg atoms b y (pul sed) electric field ionization (cf. 

e.g. laiililililiiiailiii ). 



In this paper we propose a statistical method of de- 
tecting and analyzing the physical properties of the Ry- 
dberg crystallization phenomenon and discuss its predic- 
tive power. The idea roots in the experimental procedure 
itself. A typical measurement cycle starts with the gen- 
eration of an ultracold atomic cloud and the subsequent 
excitation of a fraction of the atoms to a Rydberg state. 
Afterwards, measurements are performed, during which 
the Rydberg atoms are eventually de-excited. Then, the 
system is ready for another preparation i^^i'^^i^^" — An im- 
portant point is that the experiments are performed with 
the same number of atoms and using the same electro- 
magnetic fields in every cycle. However, since the ar- 
rangement of atoms varies between cycles, it is necessary 
to calculate statistical averages of the observables of in- 
terest. The simplest observable is the number of Rydberg 
atoms M in the cloud. The fact that for a given M, a 
regular arrangement of the Rydberg atoms on a lattice 
minimizes their interaction energy, should be visible in 
the probability distribution P{M). As we shall show be- 
low there are indeed differences between the histograms 
for crystallized and random phases. However, the pair 
correlation function turns out to possess even higher pre- 
dictive power i^i^ 

This article is organized as follows. In Sec. |TT] we shall 
introduce the model and connect its parameters to possi- 
ble experimental setups. We shall define the relevant ob- 
servables and explain the details of our numerical imple- 
mentation. The results of the calculation together with 
a thorough analysis of the arising features is contained in 
Sec, mil Section ITVl is devoted to the summary of results. 



II. THE MODEL AND SIMULATION METHOD 

We assume that in every measurement cycle the system 
consists of N atoms located at randomly chosen positions 
ri,i = 1, . . . N, uniformly distributed over the entire sys- 
tem volume.— We focus on the case of a frozen Rydberg 
gas, where the kinetic energy is negligible. Each of the 
atoms can either be excited to a Rydberg state or stay in 



its ground state. Since the electrostatic properties do not 
depend on the details of the Rydbcrg states, each atom 
can be modeled as a two-level system, and we describe 
the ensemble as a set of N randomly arranged, interact- 
ing spin-1/2 systems. Hence, the Hamiltonian reads^ 
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where ax,z denote the Pauli matrices. This can be inter- 
preted as a generalized spin-1/2 quantum Ising model. 
il is the frequency of the exciting laser and, in the spin 
language, represents a magnetic field perpendicular to 
the quantization axis, which we chose to be the z-axis. 
The detuning, i.e., the difference between the laser fre- 
quency to the resonance frequency of the Rydbcrg state, 
is denoted by A. It corresponds to a magnetic field ap- 
plied in z-direction. The third parameter, C, indicates 
the strength of an effective interaction between excited 
atoms and causes the dipole blockade explained above4^ 
We note that we are only interested in the case A > 0, 
because otherwise it is energetically not favorable to ex- 
cite atoms. This allows us to use A as a basic energy 
scale. 

An adequate modeling of the system also requires ge- 
ometrical constraints describing the trap potentials. In 
the following we shall consider different scenarios: (i) 
a ID system with open boundaries;^ (ii) a 3D system 
with open boundaries, (iii) a ID system with periodic 
boundary conditions. Options (i) and (ii) are very natu- 
ral models for realistic experiments, but make it difficult 
to extrapolate the presented numerical results towards 
realistic system sizes. Option (iii), on the other hand, is 
perfectly suitable for the calculation of correlation func- 
tions and is easily scalable to large system lengths. The 
crucial feature of our ID model is a rather large "coordi- 
nation number" , i.e., the number of atoms which interact 
significantly with any given Rydbcrg atom. This feature 
is shared by any generic 3D realization of the system up 
to some irrelevant spatial distribution parameter. That 
is why we believe that the physics in 3D is expected to 
be very similar to our ID model. 

Just as in the actual experiments we calculate averages 
over the large number of different, randomly sampled 
atom arrangements. In every cycle the effective model 
is a long-range quantum Ising model with a set of cou- 
pling constants generated by the atom positions r^. For 
every such constellation, we determine the ground state 
and evaluate the number of Rydberg atoms Af (which 
is equivalent to the magnetization in the Ising model), 
the density profile and correlation functions. The most 
difficult task is finding the ground state. We use one 
method, numerical diagonalization of the Hamiltonian 
matrix, in three different variations. For small atom 
numbers {N ss 12) the Hamiltonian can be diagonalized 



exactly. For larger numbers of atoms (TV « 30), we trun- 
cate the Hilbert space in different ways in order to speed 
up the numerical diagonalization. On the one hand, this 
truncation can be done by only keeping states in which 
the number of Rydberg atoms M remains below a cer- 
tain threshold Af*. A different approach is to only use 
the basis states with the lowest energy expectation values 
to create an effective Hamiltonian, which is then diago- 
nalized. We checked the reliability of both procedures by 
changing the respective cutoff parameter. 

To calculate the number of Rydberg atoms and the 
correlation function we proceed as follows. For a given 
random distribution of atoms, the Hamiltonian matrix is 
expressed in a basis consisting of states in which an inte- 
ger number of atoms is in the Rydberg state while the rest 
is in the ground state. The corresponding Hilbert space 
is then truncated as explained above, and the smallest 
eigenvalue and the corresponding eigenvector (the ground 
state \GS)) are obtained numerically. The number of Ry- 
dberg atoms in the ground state is now found from 
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where N is the dimension of the truncated Hilbert space, 
Mii) = E^i(«l(l + cri^^)N) the number of Rydberg 
atoms in the basis state \i) and vi — {i\GS) is the overlap 
between \i) and the ground state. One easy way of num- 
bering the basis states is to assign a "1" to a Rydberg 
atom and a "0" to a ground state atom. In this way, 
every basis state can be uniquely mapped to the binary 
representation of a number i € Nq. 

Furthermore, we calculate the pair correlation function 



.9(r) = (PRydborg(r)pRydbcrg(O)), 
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where PRydbcrg(r)dr denotes the number of Rydberg 
atoms in a volume element dr around the point r. The 
first step is to divide the interval of possible distances 
between two atoms ([0;L/2] for periodic boundary con- 
ditions) into k subintervals of equal length. Calculating 
the ground state \GS) for a given distribution of atoms 
yields the coefficients vi. Now, we consider a single pair 
of atoms and measure their distance in the current dis- 
tribution of atoms. This distance lies within one of the 
aforementioned subintervals. To the value which is al- 
ready stored for this subinterval we now add the sum of 
all |u,;p that correspond to a basis state in which both 
of the atoms of the considered pair are in the Rydberg 
state. After repeating this procedure for every possible 
pair of atoms we start over by generating a new random 
distribution. The cumulative sum of all samples treated 
this way then gives the total correlation function for a 
given set of parameters. 




III. RESULTS AND DISCUSSION 

Let us first discuss the simplest case of a noninteract- 
ing system, (7 = 0. In tfiis case, the problem is exactly 
solvable and the number of Rydberg atoms is given by 



(4) 



In the noninteracting limit, this value is independent 
of the positions r^ of the atoms. Therefore, the histogram 
P{M) becomes trivial, P{M) = 5{M - Mq). The den- 
sity distribution po = (PRydbcrg(r)), averaged over many 
realizations, is uniform, and the pair correlation function 
g{r) = Pq is constant. 

The situation changes drastically for any nonzero C. 
Figure [T] shows the data for the density distribution 
(PRydbGrg(r)) in ID and 3D systems with open bound- 
aries. While for weak interactions the Rydberg atoms 
tend to populate the boundaries, in the case of strong in- 
teractions a sizeable fraction of the atoms is redistributed 
towards the system's bulk, indicating that long-range 
order is established in the system. As we are dealing 
with a system with long-range interactions, the numeri- 
cal complexity is determined by the number of Rydberg 
atoms and not by the precise geometry of the system. 
On the other hand, going to higher dimensions at fixed 
atom number N increases the surface/bulk ratio and thus 
makes the detection of long-range order more cumber- 
some. Therefore we shall concentrate on (quasi)-lD sys- 
tems from now on. Nonetheless, we have conducted a 
number of simulations on 3D systems and found compa- 
rable results, see e.g.. Fig. [T]for the density profile of a 
3D system. We also would like to remark that an ex- 
perimentally realizable cigar-shaped confining potential 
(with a typical radius smaller than Rb) can be very well 
be approximated by the (quasi)-lD geometry considered 
here. 

Figure [5] shows a typical result for the histogram 
P{M). For r2 = 0, the number of Rydberg atoms com- 
mutes with the Hamiltonian, so the ground state for a 
given set of position {vi} has an integer expectation value 
M of Rydberg atoms. The distribution P{M), generated 
by considering all arrangements {vi}, thus becomes a se- 
ries of discrete, weighted peaks at integer values M. For 
small < r2 ^ A, these peaks broaden up, but the 
distribution remains more or less discrete. Increasing n 
leads to further broadening until eventually a continuous 
distribution is approached, see inset of Fig. [2] 

Surprisingly, in most cases the envelope of the his- 
togram can fairly well be fitted by a Gaussian, as op- 
posed to the Poissonian used, e.g., in Ref. ^^. Nonethe- 
less, the distribution function's higher order cumulants 
are not exactly zero and depend on C, which indicates 
a slight deviation from the Gaussian distribution. The 
mean fi and the variance a for a series of simulations are 
plotted in Fig. [31 Very interestingly, the fJ.{C) depen- 
dence is given by a power law with an exponent which 
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FIG. 1: Density distribution of Rydberg atoms, darker color 
indicates higher density. Left panel: ID density as a function 
of interaction strength. The system boundaries are preferred. 
Strong interactions produce a single peak in the center. Pa- 
rameters: f2/A = 0.1 and N = M* = 6 atoms (no cutoff, 
exact diagonalization) . Right panel: 2D projection of a 3D 
system with iV = M* = 6, fi/A = 4 and G = 10 in a cube 
with open boundaries. There is a low density at the center of 
the volume while it is high in the corners. 
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FIG. 2: Main figure: Histogram of the number of Rydberg 
atoms with Gaussian fit. The data points shown are the cu- 
mulative sum of all bins in the interval [n — 0.5, n + 0.5]. 
The features of the plot are explained in the text. Parame- 
ters: n/A = 0.1, C/iAL^) = 4 ■ 10"^ with iV = 10 = M* 
atoms (which corresponds to no cut-off). Inset: Histogram 



for n/A 
atoms. 



0.64, C/iAL" 



5-10"' with iV = 10 = M* 



changes abruptly from « —0.05 to ~ —0.1 at around 
C/{AL^) « 10~^, where L is the system length. We find 
this change to be a harbinger of a phase transition in the 
system. Both the mean and the width of the distribution 
function decrease as a function of C because a stronger 
repulsion increases the blockade radius Rb- This quali- 
tative behavior is independent of the value of fl. 

If one plots the mean and the variance as functions of fl 
the general behavior of the mean n{^) is qualitatively the 
same as for the noninteracting case [^(51) = //(fi, C = 0); 
cf. Eq. dlD] for n<n*, where n* = max(G/L6, A) is the 
largest energy scale in the system. For J7 > fi* the spin 
fiip term dominates the Hamiltonian ([T]) and M tends to 
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FIG. 3: Average number of Rydberg atoms n (left) and vari- 
ance a (right) in dependence on C (double- log- plots). The 
straight lines in the average number of Rydberg atoms plot 
are guides to the eye only. The parameters for all three plots 
are fl/A = 0.2 with iV = 10 = M* atoms. The resulting 
Mandel parameter is displayed in Fig. [l] 
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FIG. 5: Left: Correlation function with adjustment of the 
system length in such a way that main and secondary peaks 
coincide, C/(AL ) — 7.5 ■ 10~ . Right: Correlation function 
as a function of distance. The aliasing can be seen in the 
shoulders of the second and third maximum, C/{AL^') = 4 • 
10"®. For both correlation functions A'^ = 30 and fi/A = 25. 
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FIG. 4: Mandel Q parameter in dependence on the interaction 
strength C. The data points are calculated using the data 
from Fig. [3 




FIG. 6: Schematic illustration of the effect shown in Fig. [5] 
The dashed line is part of the correlation function in the right 
interval shifted by the length of the interval. The sum of 
dashed and solid curve in the left interval is a curve with 
main peaks that have secondary peaks on their shoulders. The 
secondary peak on the shoulder of the first peak is suppressed 
because of the blockade phenomenon. 



N/2. The only effect of C is the change in the overall 
amplitude of the curve. 

To quantify the deviation of P{M) from a Poissonian 
distribution we use the Mandel Q parameter which is 
defined as^ 






(5) 



It is plotted in Fig. [4] for the same series of simulations. 
Being negative throughout, it indicates that the P(M) 
distribution is sub-Poissonian J^i^ As realized in Ref. [2^ 
this points towards an efficient Rydberg blockade. 

Another interesting quantity is the pair correlation 
function g{r), which is defined as the probability to find 
a Rydberg atom at a distance r from another Rydberg 
atom. In order to exclude boundary effects, we now 
switch to periodic boundary conditions. To be able to 
handle computations with larger numbers of atoms we 
change the truncation procedure to one in which we con- 
sider only a fixed number of basis states. This enables us 
to approach N « 30. This set of states is chosen to be the 
one with the smallest diagonal elements in the Hamilto- 
nian matrix. We thoroughly investigated the effect of the 
truncation by considering the same system with different 
numbers of contributing states. We find that the result 
is almost independent of the number of states as long as 
it exceeds a certain threshold. All plots displayed in this 
paper meet this requirement. 



A typical correlation function is shown in Fig. [5] and is 
qualitatively comparable to the ones shown in Ref. [la |. 
On the right panel one can sec additional smaller peaks, 
which originate from the periodic boundary conditions, 
on the left side of two principal peaks. This feature, 
which is known as "aliasing" , arises in the following way: 
in the plots shown in Fig. [5] (and in any other plot of 
a correlation function in this work) our domain of defi- 
nition is the interval [0,L/2] since L/2 is the maximum 
distance between two atoms in a ID system with periodic 
boundary conditions. Now we could extrapolate this cor- 
relation function for larger distances. Since we are not 
able to resolve those the correlation function is mapped 
onto the given interval periodically. So the additional 
peaks appearing here are basically the 5th and 6th or- 
der peaks of the correlation function. Fig. [S] illustrates 
this behaviour, which is also known from the numerical 
realization of Fourier transforms with finite frequency 
cutoff4i This effect can be hidden when one uses such 
parameter values at which L and Rb are commensurate, 
see the left panel of Fig. [5l Using this control prescrip- 
tion we have calculated g{r) for a very large number of 
samples for a number of different interaction strengths, 
see Fig. [9] 

As expected, one can immediately identify the block- 
ade radius Rg as the position of the first maximum. It is 
the distance between two atoms up to which it is disad- 
vantageous to excite both of them to the Rydberg state. 
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FIG. 7: Main graph: Rb as a function of C (double-log-plot). 
The data points clearly show a power law. The parameters 
are: fl/A = 0.1 with A'' = 8 = M* atoms. The fit is discussed 
in the text. Inset: Rb as a function of Q (double-log-plot). 
The parameters are: C/{AL^) = 10"^ with N = 7 = M* 
atoms. 
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FIG. 8: Main graph: Correlation function as a function of 
distance. The dashed line represents a fit explained in the 
text. Inset: The same plot as in the main plot shown on a 
logarithmic scale. The parameters of the plot are: Q.//S. = 1, 
C/(AL«) = 2 ■ 10-5 with iV = 6 = A/* atoms. 



Obviously, it is nonzero for all interaction strengths C. 
The actual Rb{C) dependence is very interesting and is 
plotted in Fig. [7l Very naturally, Rb grows with inter- 
action strength as 



Rb « (C/AL 
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where 7 ~ 1/6 ± 1%. This was also discussed in 
Refs. 13ll3al48ll49| . Furthermore, we find the asymptotic 
behavior of the correlation function for r ^ Rb to be 
highly universal (see Fig. [8|). It is given by 



g{r) oc r 



r <^Rb 
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which is in contrast to the step-like behavior found in 
Ref. ^M,- Beyond the point r = Rb the qualitative 
shape of g{r) depends strongly on C . While for weak 



FIG. 9: Correlation function of ID samples measured in arbi- 
trary units. The curves correspond to C/AL =4-10- (dot- 
ted), C/AL** = 4-10-^ (dashed) and C / AL^ = 4-10-^ (solid) 
where all other parameters remained the same: Sl/A = 1 and 
N — 25 with 300 basis states. The curves indicate that there 
is a critical value for C/AL® between the two larger values 
given above. 



interactions, C/{AL^) < 4 • 10^^, no additional peaks 
can be seen, long-range order emerges for strong inter- 
actions C/{AL^) > IQ-^, and manifests itself as a se- 
ries of equidistant peaks. In the former case, the Ry- 
dberg atoms remain in a gas phase, whereas the latter 
situation might by described as the Rydberg crystal pro- 
posed in Refs. [2lll50l | . In an ideal crystal, the additional 
peaks would be sharp. In our calculations this cannot be 
achieved due to finite-size effects. From our simulations 
we estimate the dimensionless critical parameter for the 
transition as Ccrit/(AL^) « 5 • 10"'^. 

In fact, one can recognize this critical value for the 
phase transition already in the simple statistical param- 
eter of the P{M) histograms. Indeed, not only the mean 
value fi, but also the variance and Mandel parameter Q 
remarkably change their behavior near the critical value 
Ccrit- For instance, Q is roughly constant for C < Cent 
but starts to decrease just below Ccrit, see Fig. [4l The 
same happens to the variance. The behavior of /x is less 
prominent but still clearly detectable. Similar behavior 
can be seen in higher cumulants of P(M). However, the 
larger statistical errors might make them less useful in 
practical applications. 



IV. CONCLUSIONS 

Using exact numerical diagonalization and approxima- 
tive descendants of this method we have investigated the 
Rydberg crystallization phenomenon in ultracold gases. 
We have estimated the critical interaction strength by 
two different techniques. Both pair correlation function 
and simpler statistical data show signatures of this phase 
transition. The big advantage of using mean and vari- 
ance of the histogram for the number of excited atoms 
measured in a long series of identical experiments is its 



good experimental accessibility. In this way we have de- 
veloped a purely statistical method of detecting Rydberg 
crystallization. We hope that this procedure can soon be 
implemented in state-of-the-art experiments in order to 
unambiguously identify the Rydberg crystallization phe- 
nomenon without complicated spectroscopic techniques. 
Furthermore, the presented details of the pair correla- 
tion function might be useful for the continuous experi- 
mental efforts in spatial imaging of Rydberg aggregates 
in the spirit of Refs. |36ll4ll |. This technique is not only 
able to yield very precise value of the blockade radius, but 
has also proven to be a reliable source for the estimation 



of the critical parameters for Rydberg crystallization. 
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